ENERGY  RELEASE  FROM  CONDENSED  PHASE 
MATERIALS  AND  HETEROGENEOUS  REACTIVE 
FLOW  MODELING 
FINAL  REPORT 
SAI-83-831-WA 
Prepared  by 

Kazhikathra  Kailasanath  and  Ellis  Hyman 


ENERGY  RELEASE  FROM  CONDENSED  PHASE  MATERIALS 
AND  HETEROGENEOUS  REACTIVE  FLOW  MODELING 

SAI -  83-  831- KA 

Final  Report 

Submitted  to: 

Laboratory  for  Computational  Physics 
Naval  Research  Laboratory 
Washington,  D.C.  20375 

Prepared  Under 

Contract  No.  N00014- 81-C- 2418 

Prepared  by 

Kazhikathra  Kailasanath 
and 

Ellis  Hyman 


SCIENCE  APPLICATIONS,  INC. 

1710  Goodridge  Drive,  P.O.  Box  1303 
McLean,  Virginia  22102 
(703)  734-5840 


SECURITY  CLASSIFICATION  OF  This  FAGS  <TWiw  Bate  gntered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

'l.  NUMRER  2.  OOVT  ACCSSSION  NO. 

/77 

1.  RECIFiENT'S  CATALOG  NUMRER 

4.  TITLE  fund  Subtitle) 

ENERGY  RELEASE  FROM  CONDENSED  PHASE 
MATERIALS  AND  HETEROGENEOUS  REACTIVE 

FLOW  MODELING 

S.  TYFE  OF  REFORT  S  FERtOD  COVERED 

Final  Report 

18  Jun  81  -  17  Jun  82 

7.  AUTHOR?  •> 

Kazhikathra  Kailasanath  and  Ellis  Hyman 

S.  CONTRACT  OR  GRANT  NUMBER?*; 

N00014-81-C-2418 

*■  FERFONMINO  ORGANIZATION  name  ano  aooress 

Science  Applications,  Inc. 

1710  Goodridge  Drive 

McLean,  VA  22102 

10.  frooram  ELEMENT.  FROjLCT.  task 
AREA  S  WORK  UNIT  NUMBERS 

II.  CONTROLLING  OFFICE  NAME  AND  AOORESS 

Naval  Research  Laboratory 

4555  Overlook  Avenue,  S.W. 

Washington,  D.C.  20375 

11.  REFORT  OATE  ’ 

July  1982 

M.  MONlVoRiNG  AGENCY  NAME  *  ADDRESS?//  dt  Hu  tent  hum  Controlltnd  Olll eu) 

IS.  SECURITY  CLASS,  (ul  Hi la  re pert) 

UNCLASSIFIED 

IS.  OISTRIRVlTlON  STATEMENT  (ol  IMa  RUbUrt) 

DISTRIBUTION  STATEMENT  A 

Approved  Joe  public  release 
Distribution  Unlimited 

_ L 

17.  DIST  R1SUTION  STATEMENT  ft  thm  tract  in  Bt—k  10,  II  BltHrmi  Ns  / 

IS.  SUFFLEMENTARY  NOTES 

IS.  KEY  WOROS  (Camthme  un  imm  NR  II  Mltiwry  and  Identity  Uy  black  mbNtJ 

Energy  Release  Reactive  Flow 

Flames  Droplet  Combustion 

Detonations 

20.  ABSTRACT  (Cumlmn  un  tuaueee  aid*  If  nimm)  and  Identity  hr  MM  number) 

This  report  describes  the  results  from  investigations  on 
(1)  the  burning  velocity  of  hydrogen  in  air,  (2)  the  direct 
initiation  of  detonations  in  hydrogen-air  mixtures  and  (3)  the 
flow  in  and  about  fuel  droplets. 

A  time  dependent  numerical  model  has  been  used  to  simulate 
a  flame  propagating  in  a  stoichiometric  hydrogen-air  mixture. 

(■Continues) 

00  I  jmTtS  1473  COITION  OF  1  NOV  •>  IS  OMOLETE 

S/N  01 03- LF -01 4.4401 


SECURITY  CLASSIFICATION  OF  THIS  FAOE  (Km  SSa  Rm« 


SECURITY  CLASSIFICATION  OF  THIS  face  (Ska  Omm  bMn« 


20.  Abstract  (Continued) 

The  burning  velocity  of  hydrogen  in  air  was  estimated  from  the 
numerical  simulations  and  was  found  to  be  within  the  observed 
range  of  experimental  values.  The  effect  of  inaccuracies  in 
the  estimation  of  the  molecular  diffusion  coefficients  and  in 
the  initial  conditions  have  also  been  discussed. 

Numerical  simulations  of  shock  tube  experiments  have  been 
used  to  show  that  a  minimum  energy,  a  minimum  power  and  some 
relation  between  them  is  required  to  characterize  the  direct 
initiation  of  detonations.  These  simulations  also  show  that 
under  certain  conditions ,  the  induction  time  of  the  gaseous 
mixture  studied  is  very  sensitive  to  perturbations  in  pressure 
and  temperature.  Simulations  have  been  carried  out  to  evaluate 
the  effects  of  sound  wave  perturbations  which  in  turn  result  in 
temperature  and  pressure  fluctuations.  The  results  of  these 
sound  wave  studies  show  how  the  ignition  time  is  affected  by 
sound  waves  of  a  given  amplitude  and  frequency. 

The  initial  results  from  numerical  simulations  of  the 
inviscid,  incompressible  flow  in  and  about  a  non-reacting  fuel 
droplet  have  also  been  presented. 
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GENERAL  DISCUSSION 


During  this  contract  year  we  have  performed  a 
number  of  computer  simulations  to  investigate  (1)  the 
burning  velocity  of  hydrogen  in  air,  (2)  the  direct  initia¬ 
tion  of  detonations  in  hydrogen-air  mixtures,  and 
(3)  compressible  flow  around  fuel  droplets.  The  efforts 
and  accomplishments  in  each  of  these  areas  will  be  described 
in  some  detail  in  the  following  paragraphs. 

We  estimated  the  burning  velocity  of  hydrogen  in 
air  using  a  time  dependent  model  for  a  flame  burning  in 
a  hydrogen-oxygen-nitrogen  mixture,  and  compared  our 
results  with  experimental  data.  In  the  past,1  SAI ,  in 
collaboration  with  NRL  has  used  the  NRL  ID  flame  code  to 
study  the  ignition  and  quenching  of  flames  in  pre-mixed 
gases.  The  flame  code  has  now  been  extended  to  study  the 
propagation  of  laminar  flames  in  pre-mixed  gases.  It  was 
found  that  the  time-dependent  flame  model  is  capable  of 
predicting  both  the  transient  and  steady  state  properties 
of  laminar  flames  in  hydrogen-oxygen-nitrogen  mixtures. 

The  code  was  used  to  estimate  the  burning  velocity  of 

hydrogen  in  air  and  this  has  been  compared  with  experimental  data . 

The  burning  velocity  of  in  air  at  1  atm  pressure 
and  298°K  predicted  by  the  model  was  at  the  lower  end  of 
the  observed  range  of  experimental  values.  We  note,  how¬ 
ever,  that  there  are  several  factors  which  would  tend  to 


increase  this  estimated  value.  First,  the  calculated 
burning  velocity  would  be  larger  if  the  velocity  of  the 
unburned  gas  used  in  this  calculation  was  closer  to  the 
flame  front  instead  of  at  the  open  boundary.  It  is  not 
clear  where  it  should  be  evaluated.  Second,  thermal  dif¬ 
fusion  effects  have  not  been  considered  in  the  model. 

Finally,  the  extremely  important  hydrogen  diffusion  coef¬ 
ficients  are  not  really  known  at  high  temperatures.  If 
they  are  higher  than  the  ones  used  here,  the  flame  velocity 
would  increase.  These  observations  lead  us  to  believe 
that  the  higher  observed  values  of  the  burning  velocity 
are  probably  more  accurate. 

This  study  has  also  brought  out  the  importance 
of  initial  conditions.  As  in  steady-state  calculations, 
a  better  initial  guess  leads  to  faster  convergence  to  a 
steady  propagating  flame.  Perhaps  the  best  way  to  initialize 
such  a  problem  is  to  perform  a  two  step  process  in  which 
we  first  determine  the  minimum  ignition  energy  for  a  given 
radius  and  time  of  energy  deposition  (as  reported  in 
Ref.  1),  and  then  deposit  this  energy  and  evaluate  the 
properties  of  the  propagating  flame  front.  However, 
considering  that  the  experimental  velocities  vary  by 
~25l,  a  good  initial  guess  for  the  temperatures  and 
species  concentrations  certainly  gives  an  adequate  esti¬ 
mate.  Finally,  we  note  that  the  criterion  on  which  to 
judge  the  initial  guess  values  is  to  look  at  how  flat  the 
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temperature  profile  is  behind  the  flame  front.  A  bad 
guess  produces  some  structure  and  oscillations  which  only 
slowlv  disappear.  This  work  has  been  presented  at  the 
' ' G AMM -Workshop  on  Numerical  Methods  in  Laminar  Flame 
Propagation"  and  appears  in  this  report  as  Appendix  A 
entitled,  "Time-Dependent  Simulation  of  Flames  in  Hydrogen- 
Oxygen-Nitrogen  Mixtures."  It  will  appear  in  the  Proceedings 
of  the  G AMM- Workshop  on  Numerical  Methods  in  Laminar  Flame 
Propagation,  Volume  6,  in  the  series  Notes  on  Numerical 
Fluid  Dynamics, 

To  investigate  the  direct  initiation  of  detonations 
in  hydrogen-air  mixtures  we  performed  numerical  simulations 
of  shock  tube  experiments.  Whereas  these  simulations 
qualitatively  agree  with  earlier  experimental  observations*'2,3^ 

in  that  a  minimum  power,  a  minimum  energy  and  some  relation 
between  them  is  required  to  characterize  the  direct  initia¬ 
tion  of  detonations,  the  experiments  and  calculations  do 
not  agree  well  quantitatively.  During  the  course  of  this 
study,  it  was  found  that  the  disagreement  is  primarily 
because  the  induction  times  of  the  gaseous  mixtures,  at 
the  temperatures  and  pressures  expected  behind  the 
incident  shocks,  are  very  sensitive  to  perturbations  in 
pressure  and  temperature. 

Simulations  were  carried  out  to  evaluate  the 
effects  of  sound  wave  perturbations  which  in  turn  result 
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in  temperature  and  pressure  fluctuations .  The  results  of 
these  sound  wave  studies  showed  how  the  ignition  time  is 
affected  by  sound  waves  of  a  given  amplitude  and  frequency. 
It  was  observed  that  low  frequency  perturbations  (~480  Hz 
or  half  wavelength  L  ~  75  cm)  of  large  enough  amplitude 
(~200  m/s)  can  reduce  the  ignition  time  by  an  order  of 
magnitude.  It  was  found  that  low  frequency  perturbations 
were  more  likely  than  high  frequency  perturbations  to  cause 
the  observed  discrepancies  between  the  experiments  and 
calculations . 

The  presence  of  temperature  and  pressure  variations 
in  the  experiments  can  be  due  to  factors  such  as  the  non¬ 
ideal  break-up  of  the  diaphragm,  boundary  layers  causing 
non-uniformities  in  the  walls  of  the  tube,  the  sudden 
change  in  shape  and  area  of  cross-section  between  the 
driver  and  driven  sections,  or  heat  loss  to  the  walls. 
However,  perturbations  required  to  substantially  alter  the 
ignition  times  and  cause  discrepancies  such  as  observed 
here  are  larger  than  what  might  normally  be  present  in 
very  well  controlled  shock  tube  experiments.  This  needs 
to  be  investigated  before  one  can  confidently  use  shock 
tube  studies  to  determine  the  critical  conditions  necessary 
to  characterize  the  direct  initiation  of  detonations. 

A  detailed  report  of  this  study  is  enclosed  as 
Appendix  B.  This  report  which  is  to  be  submitted  for 


publication,  is  entitled,  "Direct  Initiation  of  Gaseous 
Detonations  in  Shock  Tubes:  Numerical  Simulations  of 
Experiments . " 

SA1  in  cooperation  with  NRL  has  initiated  the 
development  of  a  reactive  flow  Lagrangian  triangular  grid 
hydrodynamic  code  to  treat  the  problem  of  compressible 
flows  past  deformable  droplets.  This  development  began 
by  modifying  SPLISH,  a  Lagrangian  triangular  code  for 
incompressible  fluid  flow'  in  two  dimensional  Cartesian 
geometry.  This  code  had  previously  been  used  to  study 
Rayleigh-Taylor  instabilities  and  flows  past  solid  objects. 
Our  initial  modification  was  to  convert  the  solid  object 
to  a  denser  fluid  by  tessalating  the  object  and  then 
allowing  the  code  to  treat  the  object  as  a  fluid  of  differing 
density.  The  triangular  grid  allows  us  to  place  sides  to 
triangles  to  approximate  the  droplet  fluid  interface  while 
the  Lagrangian  nature  allows  us  to  track  these  interfaces 
while  the  droplet  deforms  in  a  background  flow. 

A  typical  calculation  of  this  type  is  shown  in 
figures  1  and  2.  For  this  calculation  a  flow  from  left  to 
right  is  initialized  by  a  pressure  pulse  at  the  left 
boundary.  For  all  later  times  the  right  and  left  sides 
are  subjected  to  periodic  boundary  conditions.  The  droplet 
density  to  background  fluid  density  ratio  for  this  run  is  2. 

Figure  1  shows  the  triangular  mesh  at  several 
times  during  the  calculation.  Early  in  the  calculation  we 
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see  fluid  flows  typical  of  flows  past  solid  objects.  A 
recirculation  cone  forms  in  the  wake  of  the  droplet.  The 
droplet  eventually  compresses  in  the  direction  of  the  flow 
and  is  forced  out  normal  to  the  flow.  The  shear  flow  then 

pulls  the  top  and  bottom  of  the  droplet  around  the  recircula¬ 

tion  cone.  As  the  droplet  thins  around  this  recirculation 
cone  it  begins  to  break  into  several  smaller  pieces.  Fig¬ 
ure  2  shows  the  pathlines  of  the  vertices  for  this  same  cal¬ 
culation.  Subtracting  the  mean  flow  would  show  a  large 
recirculating  double  vortex.  Calculations  with  larger  mass 
ratios  show  similar  behavior  except  that  the  larger  inertia 
of  the  droplet  makes  the  shearing  around  the  recirculation 
cone  more  difficult. 

In  order  to  get  closer  to  more  realistic  flows  we 

have  produced  and  are  currently  testing  an  algorithm  to 

simulate  surface  tension  effects.  Now,  the  finite  difference 
formulae  used  by  SPLISH  for  the  hydrodynamics  obey  the 
conservation  laws  associated  with  the  differential  equations. 
We  have  formulated  the  surface  tension  forces  as  the  gradient 
of  a  "potential"  which  acts  only  at  the  interface.  In  this 
way  we  can  ensure  that  the  finite  difference  algorithms 
will  still  be  conservative. 

The  surface  tension  potential  is  given  by  Laplace's 
formula  for  the  jump  in  pressure  across  an  interface 
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where  P^  is  the  pressure  just  inside  the  droplet  at  the 
interface,  PQ  is  the  pressure  just  outside  the  droplet  at 
the  interface,  a  is  the  surface  tension  coefficient  asso¬ 
ciated  with  the  two  media  which  define  the  interface  and  R 
is  the  radius  of  curvature  of  the  cylindrical  droplet.  The 
radius  of  curvature  is  positive  at  points  on  the  interface 
where  the  droplet  surface  is  convex  and  negative  when  the 
surface  is  concave.  These  pressure  jumps  are  included  in 
the  Poisson  equation  for  the  pressure.  The  average  pres¬ 
sure,  (Pi+PQ)/2,  is  computed  at  an  interface  vertex.  From 
the  average  pressure  and  the  pressure  jump  we  can  compute 
a  pressure  gradient  centered  on  triangles  for  inclusion  in 
the  momentum  equation. 

The  radius  of  curvature  is  computed  from  a 
parametric  cubic  spline  interpolant  to  the  interface  ver¬ 
tices.  The  parameter  is  a  pseudo-arc  length,  s,  defined 
by  the  Spline  Knots 


where  r^,  i  =  1,  ...,  N  are  the  N  vertices  which  define 
the  interface  with  r^  =  r^ .  The  curvature  is  then  defined  by 
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K  =  | 1/R| 


l(r'x?  )  •  e„| 


where  the  prime  indicates  differentiation  with  respect  to 
the  parameter  s.  The  sign  of  R  at  an  interace,  r ^ ,  is 
given  by  the  sign  of  x  (r^-r^_j)  •  e_  where  e_ 

is  the  unit  vector  in  the  z  direction  (the  direction  normal  to 
the  interface) . 

We  have  tested  this  algorithm  by  computing  the 
period  of  a  droplet  undergoing  oscillations  due  to  surface 
tension  forces.  In  typical  runs  we  found  periods  that  were 
approximately  15%  larger  than  those  predicted  by  linear 
theory.  Similar  errors  were  found  for  Rayleigh-Tavlor 
oscillations  and  decreased  with  more  refined  meshes. 

These  results  were  reported  on  at  the  Eastern 
States  Combustion  Meeting  in  Pittsburgh,  PA  in  October  1981. 
An  abstract  of  that  paper  follows: 


Calculations  of  Flow  in  and  About  Fuel  Droplets, 

M.J.  Fritts,  D.E.  Fyfe*,  and  E.S.  Oran,  Naval  Research 
Laboratory.  Spray  combustion  is  complicated  by  a  parti¬ 
cularly  rich  variety  of  physical  processes  which  take 
place  over  a  wide  range  of  space  and  time  scales.  One 
of  the  important  aspects  of  the  problem  is  the  individual 
and  collective  behavior  of  burning  droplets.  However, 
including  the  detailed  interactions  between  droplets 
and  the  external  flow  field  into  combustion  models 
greatly  increases  the  computational  complexity.  To 
reduce  the  cost  of  these  calculations  approximations 
must  be  made  that  are  particularly  severe.  The 
important  assumptions 1 » 2  fall  into  two  major  areas: 
those  concerned  with  droplet  behavior  and  those 
concerned  with  the  details  of  the  external  flow  field. 
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The  shape  of  the  droplets  in  spray  combustion  models 
are  typically  assumed  to  be  spherical,  and  empirical 
expressions  are  used  to  account  for  drag  and  convection. 
Since  the  droplets  are  known  to  deform,  an  equivalent 
sphere  is  used  to  approximate  the  assumed  deforma¬ 
tion.  The  effects  of  droplet  breakup  are  included  by  using 
estimated  breakup  times  and  drop  sizes  after  breakup. 

Similar  assumptions  are  made  concerning  the  flow 
about  the  droplets.  Typically  the  flow  is  considered 
to  be  quasisteady.  Thus  the  alteration  of  the  flow 
field  due  to  deformations  and  wake  effects  are 
neglected,  as  is  the  interaction  of  the  flow  field 
back  onto  the  droplet.  Lack  of  knowledge  about  the 
flow  field  extends  to  larger  scales  as  well:  in  most 
models  the  droplet  concentration  is  assumed  to  be 
dilute  since  little  is  known  about  droplet-droplet 
or  droplet-wake  interactions.  Furthermore,  the 
influence  of  transport  in  the  external  flow  is  based 
solely  on  ambient  properties  because  coherent  or 
turbulent  features  pose  too  large  a  complication. 

Once  these  basic  assumptions  about  the  droplet  and 
exterior  flow  have  been  invoked,  other  restrictions 
result  because  of  the  lack  of  information  about  the 
detailed  flow  fields  both  internal  and  external  to 
the  droplet.  For  example,  transport  properties 
within  the  drop  can  be  treated  in  limiting  cases 
only,  chemical  reactions  are  neglected  in  the  liquid 
phase  and  in  the  exterior  flow  field,  the  properties 
of  the  flow  field  are  assumed  to  be  constant  at  each 
instant  in  time  with  only  concentration  diffusion 
considered,  and  pressure  is  assumed  constant  and 
equal  to  the  local  ambient  pressure. 

In  this  presentation  we  will  describe  the  results  of 
a  new  research  program  aimed  at  providing  some  under¬ 
standing  of  the  behavior  of  the  internal  and  external 
flow  fields  in  droplet  combustion  problems.  One  goal 
of  the  program  is  to  help  clarify  the  utility  and 
appropriateness  of  some  of  the  assumptions  listed 
above.  The  program  is  centered  on  a  conservative, 
Lagrangian  description  of  the  combustion  process. 

This  is  for  several  reasons. 3  First,  the  hydro¬ 
dynamics  becomes  simpler  in  that  material  boundaries 
evolve  naturally  and  are  followed  without  further 
approximations.  If  the  method  conserves  vorticity 
as  well,  these  material  boundaries  will  be  convected 
properly  by  any  coherent  vortex  patterns  which  evolve 
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in  the  external  flow  field.  Second,  calculations 
of  the  combustion  processes  themselves  become  simpler, 
since  local  fluid  parcels  are  followed  self-consistently 
through  their  changing  environment.  The  size  of  heat 
release  in  the  flow  field  is  easily  determined  locally, 
as  are  its  effects  on  neighboring  fluid  elements. 

Third,  several  numerical  advantages  are  incurred,  such 
as  elimination  of  the  need  to  difference  advective 
terms  and  the  elimination  of  numerical  diffusion  result¬ 
ing  from  lack  of  knowledge  of  interface  positions. 

The  Lagrangian  treatment  used  in  the  program  is  based 
on  a  finite  difference  treatment  using  an  irregular, 
dynamically  restructuring  triangular  grid.^  A  tri¬ 
angular  grid  avoids  problems  of  mesh  tangling  encount¬ 
ered  in  Lagrangian  implementations  using  a  quadri¬ 
lateral  mesh.  Individual  mesh  points  are  reconnected 
to  allow  for  the  migration  of  fluid  elements  in  the 
flow  field.  Since  the  number  of  grid  lines  meeting 
at  a  vertex  is  variable,  the  resolution  can  be 
altered  where  needed  (e.g.,  around  the  droplet  in 
the  combustion  zone)  without  affecting  the  resolution 
in  other  areas  of  the  computation.  Integral  equations 
are  used  to  specify  the  values  of  physical  variables 
during  grid  restructuring  to  ensure  that  the  evaluation 
of  the  fluid  flows  are  fully  conservative. 

In  this  presentation  we  describe  calculations  of 
inviscid,  incompressible  fluid  flow  about  non-reacting 
droplets  for  a  variety  of  flow  speeds  and  droplet 
densities.  The  evolution  of  the  flow  field  external 
and  internal  to  the  droplet  will  be  shown.  Preliminary 
results  which  include  the  effects  of  surface  tension 
will  also  be  shown.  The  algorithms  for  a  particularly 
simple  extension  to  compressible  flows  will  be  discussed 
for  the  case  in  which  heat  release  occurs  on  a  time 
scale  fast  compared  to  the  fluid  flow.  When  the  time 
scales  are  comparable,  or  if  the  effects  of  shocks  are 
important,  other  algorithms  must  be  used,  and  alternate 
formulations  will  be  presented. 

Future  work  on  this  project  will  incorporate  the 
compressibility  algorithms  into  the  code  together  with 
the  effects  of  viscosity.  Algorithms  for  calculating 
reactions,  vaporization,  thermal  conduction  and  molecular 
diffusion  will  be  taken  from  other  reactive  flow  codes. ^ 
The  addition  of  each  of  these  processes  into  the  code 
will  be  accompanied  by  benchmarking  on  pertinent 
physical  problems.  An  r-z  version  of  the  code  is  being 
developed  in  parallel  with  the  Cartesian  version,  and 
initial  tests  of  this  code  using  circular  Couette  flow 
and  Taylor  vortex  flow  will  be  presented  as  well. 
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Numerical  Calculations  of  Flow  in  and  around  a  Droplet. 

D.E.  Fyfe*,  M.J.  Fritts,  and  E.S.  Oran,  Naval  Research 
Laboratory.  Complicated  flow  fields  arise  in  and 
around  burning  fuel  droplets.  Here  we  present  the 
initial  results  of  numerical  simulations  aimed  at 
providing  a  description  of  these  flows.  These  initial 
calculations  show  the  results  of  inviscid  incompres¬ 
sible  fluid  flow  in  and  about  a  non-reacting  fuel 
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droplet.  The  simulations  were  performed  using  a 
Lagrangian  algorithm  with  a  dynamically  restructuring 
triangular  grid.  These  algorithms  allow  variable 
sized  cells,  addition  and  deletion  of  cells,  and 
continuous  calculation  without  diffusive  rezoning. 

We  present  the  evolution  of  the  flow  fields  for 
several  flow  speeds  and  droplet  densities.  Pre¬ 
liminary  results  incorporating  the  effects  of  surface 
tension  will  also  be  shown.  More  realistic  calcula¬ 
tions  in  the  future  will  have  to  include  viscosity, 
thermal  effects,  vaporization,  molecular  diffusion, 
chemical  reactions  and  compressibility.  Algorithms 
for  the  extension  to  compressible  flows  when  heat 
release  occurs  on  a  time  scale  fast  compared  to  the 
fluid  flow  will  be  discussed. 
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FIGURE  CAPTIONS 


Figure  1 


Figure  2 


Calculation  of  droplet  deformation  in  an 
external  flow. 

Vertex  pathlines  showing  fluid  flow  inside 
and  outside  deforming  droplet. 
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Appendix  A 

TIME -DEPENDENT  SIMULATION  OF  FLAMES 
IN  HYDROGEN -OXYGEN-NITROGEN  MIXTURES 


Time -Dependent  Simulation  of  Flames  in  Hydrogen-Oxy gen- 

Nitrogen  fixtures 

K.  Kailasanath*,  E.S.  Oran,  J.P.  Boris,  and  T.R.  Young 
Laboratory  for  Computational  Physics 
Naval  Research  Laboratory 
Washington,  D.C.,  20375 

1.  INTRODUCTION 

This  paper  describes  a  one-dimensional,  time-dependent, 
Lagrangian  model  developed  to  study  the  initiation,  propaga¬ 
tion  and  quenching  of  laminar  flames.  The  model  Incorporates 
a  number  of  new  approaches  and  algorithms  which  have  now  been 
tested  by  comparisons  to  less  complex  or  analytic  solutions 
and  by  comparisons  to  experimental  data.  These  new  elements 
include:  ADINC  [1]  an  implicit,  Lagrangian  method  for  solving 

the  convective  parts  of  the  conservation  equations;  DFLUX 
[2,3]  a  variable  accuracy  algorithm  for  determining  diffusion 
fluxes  without  having  to  invert  matrices;  SPLIT  and  MERGE, 
routines  for  dividing  or  merging  computational  cells  as 
specified  by  external  criteria;  VSAIM,  a  vectorized  version  of 
the  ordinary  differential  equation  solver,  CHEMEQ  [*1,5];  and  a 
new  method  for  treating  an  open  boundary  In  an  implicit, 
Lagrangian  calculation.  An  asymptotic  coupling  method,  used 
in  conjunction  with  timestep  splitting  to  couple  the  various 
processes,  allows  the  use  of  entirely  different  algorithms  for 
the  physical  processes  represented  by  different  mathematical 
forms. 

The  model  has  been  used  for  a  variety  of  flame  studies  of 
hydrogen-oxygen -nitrogen  mixtures.  These  include  calculations 
of  minimum  ignition  energies,  flammability  limits,  quench 
volumes,  and  burning  velocities.  The  chemical  rate  scheme  has 
now  been  tested  extensively,  as  have  the  thermal  and  molecular 
diffusion  coefficients.  Thus  we  expect  the  model  to  calculate 
correctly  the  time-dependent  behavior  Implied  by  the  Initial 
and  boundary  conditions  supplied. 

*Currcntly  with  Science  Application,  Inc.  McLean,  VA 


After  a  summary  is  given  of  the  important  numerical 
features,  several  calculations  are  presented.  The  first  is  a 
flame  initiation  and  minimum  ignition  energy  study  of  a 
mixture  of  H2 :02 :N2/2: 1 : 10  initially  at  298  K  and  1  atm. 

These  are  inherently  time-dependent  problems  and  take  advan¬ 
tage  of  this  property  of  the  model.  The  second  problem  dis¬ 
cussed  is  a  calculation  of  the  burning  velocity  of  the  mixture 
H2  :C  :K2/2: 1:4,  again  at  298  K  and  1  atm.  The  burning 
velocity  is  a  quantity  used  to  describe  a  steady  state 
property  and  it  is  generally  calculated  by  steady  state 
methods.  We  show  below  that  the  model  described  also  does 
well  on  this  type  of  calculation. 

2.  NUMERICAL  MODEL 

2. 1  Basic  Equations  ‘  ' 

We  solve  the  time-dependent  equations  for  conservation  of 
total  mass  density  p,  momentum  pv,  and.  energy  E  as  well  as  the 
individual  species  number  densities  {nj}.  These  may  be 
written  as  [2,63: 

| £«-I.pv  (1) 

an, 

aT  "  -I-njX j  +  pj  "  LjnJ  (2) 

rn 

■gfc—  -V.(pvv)  -  v  P  +  v.  nm[Vv  +  ( Vv)  3  (3) 

||  =  -v.Ev  -  v. Pv  -  V.Q  (A) 


where  the  heat  flux  vector,  Q,  is  defined  as 

T 

2  -  -V?  +  $"jh  A  *  V  j^kSwJS^  A  -  V 


The  quantity  v  is  the  fluid  velocity,  the  {Vj }  are  the  dif¬ 
fusion  velocities,  and  { P j }  and  {L j }  refer  to  the  chemical 
production  and  loss  processes  for  the  individual  species  J. 
The  quanties  nm  and  Xm  are  the  mixture  viscosities  and 
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thermal  conductivities  respectively.  The  superscript  "T"  in 
the  last  term  of  Eq.  (3)  indicates  that  the  transpose  is 
taken.  The  quantities  {h i }  are  the  temperature  dependent 
enthalpies  for  each  species,  and  the  (Djk }  ard  } 

are  the  sets  of  binary  and  thermal  diffusion  coefficients, 
respectively. 

We  also  assume  that  the  mixture  consists  of  ideal  gases 
so  that  the  pressure,  P,  may  be  written  as: 

P  =  NkgT  *  (6) 

where  N  is  the  total  number  density,  kg  is  Boltzmann's 
constant  and  T  is  the  temperature.  The  model  presented  in 
this  report,  however,  is  not  restricted  to  ideal  gases  and  in 
fact  any  equation  of  state  may  be  used. 

The  diffusion  equations  may  be  written  as 


where  the  source  terms  Sj  are  defined  as 


.  .  T(is,  (ii  Hi,  =?  .  HjHfc  S 

V  -(i!  )  ~  (p  n  }  p "  £  (pk  ~  p3  ) 


VT 

T 


(8) 


The  diffusion  velocities  are  also  subject  to  the  constraint: 

|  ’/j  - 0  (9> 


2.2  Convective  Transport 

The  convective  transport  terms  in  Eqs.  (1-*1)  are  solved 
by  the  algorithm  ADINC  Cl]  •  ADINC,  designed  for  either 
Adiabatic  or  incompressible  flows,  is  an  implicit  and 
Lagrangian  algorithm.  Since  it  communicates  compression  and 
expansion  across  the  system  implicitly,  it  overcomes  the 
Courant  time-step  limit.  Since  it  is  Lagrangian,  it  can  main¬ 
tain  steep  gradients  computationally  for  a  long  period  of 
time.  This  is  important  in  flame  calculations  where  the  dif¬ 
fusive  transport  of  material  and  energy  can  govern  the  system 
evolution  and  therefore  must  be  calculated  accurately. 

ADINC  solves  the  following  equations  for  mass  and  momen¬ 
tum  transport  in  one  dimension: 


(10) 


dv 
p  dt 


(ID 


=  -  VP 

The  energy  evolution  equation  is  eliminated  by  using  an 
adiabatic  equation  of  state: 

p(P,S)  «  p  +  (P/S)l'Y  (12) 

v 

This  equation  of  state  with  pc  =  0  is  correct  for  adiabatic 
compression  and  expansion  of  an  ideal  gas.  The  entropy,  S,  is 
assume",  constant  throughout  the  numerical  integration.  Non- 
adiabatic  processes,  such  as  external  heating,  thermal  conduc¬ 
tion  and  chemical  energy  release,  are  addeu  to  Eqs.  (10-11) 
using  the  timestep  splitting  methods  described  below. 

Given  an  aproximation  to  the  pressure,  ADINC  calculates 
the  fluid  density  using  Eq.  (12).  This  equation-of -state 
density  is  compared  to  the  density  derived  from  the  fluid 
dynamics  through  Eq.  (10).  The  difference  is  iterated  to  zero 
using  a  quadratically  convergent  implicit  solution  of  Eq.  (11) 
which  then  gives  an  improved  approximation  to  the  pressure. 
During  this  iteration  that  analytic  derivative  dA/dP  is  used 
where  A  is  th .  volume  of  a  computational  cell.  Thus 

I5F  '  -h?  (P/S)W*  <«> 

for  the  particular  equation  of  state  (12).  ADINC  also  assumes 
that  pressure  ar.d  density  are  constant  within  each  individual 
finite-difference  cell  and  that  the  physics  is  evolving  slowly 
enough  for  full  communication  across  that  cell  to  have  occur¬ 
red  in  a  timestep. 

ADINC  has  been  used  extensively  for  solving  a  wide  vari¬ 
ety  of  problems.  Some  of  these  and  a  number  of  tests  of  the 
algorithm  have  been  documented  by  Boris  [1], 


2.3  Chemical  Kinetics  Calculations 

The  coupled,  nonlinear,  ordinary  differential  equations 

which  describe  the  chemical  interactions  are  taken  from  that 

part  of  Eq.  (2),  which  represents  the  production  and  loss  of 

reacting  species: 

3n  i 

■at  “  pj  "  Lj  nJ  *  J  "  i*  *  .  ’  (1/|) 


where  M  is  the  total  number  of  species  present.  The  func¬ 
tional  dependencies  of  the  terms 

L  j  —  L  j  {n^  (t)}  >  k  -  1 ,  • .  . . ,  M  (15) 

emphasize  the  strong  coupling  between  the  various  species. 

The  system  represented  by  Eq.  (14)  may  be  stiff  when  there  are 
large  differences  in  the  time  constants  associated  with  dif¬ 
ferent  chemical  reations.  Stiffness  may  occur  for  different 
species,  in  different  locations,  at  different  times  or  simul¬ 
taneously  throughout  the  course  of  an  integration.  Because  of 
this,  when  there  are  a  large  number  of  reactions,  the  solution 
of  the  chemical  kinetics  equations  is  usually  the  most  expen¬ 
sive  part  of  a  detailed  reactive  flow  calculation.  Further¬ 
more,  the  computational  cost  increases  with  the  number  of 
species  and  the  dimensionality  of  the  problem.  Therefore  a 
method  which  is  efficient,  accurate,  conservative,  stable  and 
which  does  not  require  storage  of  large  quantities  of  data 
from  one  timeste?  to  another  is  required.  Such  a  method  is 
VS AIM,  which  is  a  fully  vectorized  version  of  the  selected 
asymptotic  integration  method  employed  in  CHEMEQ  [4,5). 

In  this  method,  the  stiff  equations  are  identified  and  solved 
using  a  very  stable  asymptotic  method.  The  remaining 
equations  are  solved  using  a  standard  classical  method. 

Once  the  species  number  densities  are  knotv'n,  the  temper¬ 
ature  can  be  evaluated  by  iterating  on  the  equation: 

Vli,!-;v"p;vr'  (16) 

J  J  J 

where  E^  is  the  total  energy  per  unit  volume,  hQj  is  the 
set  of  heats  of  formation  of  the  species  J  and  hj  is  the 
sensible  enthalpy  of  species  j.  A  Ncwton-Raphson  iteration 
technique  is  used  and  it  usually  converges  in  two  or  three 
iterations  per  timestep. 


2. 1|  Diffusive  Transport 


The  diffusive  transport  processes  considered  in  this 
model  are  molecular  diffusion  and  thermal  conduction.  These 
are  the  parts  of  Eqs.  (1-4)  which  are  represented  by  the 
expressions : 


TT5  "  -  i-n  jSj 

(17) 

’  !-E-»2T  + 

(18) 

The  above  equations  are  conservatively  differenced  and  solved 
explicitly.  The  effects  of  viscosity  and  thermal  diffusion 
have  not  been  considered  for  the  results  discussed  in  this 
paper. 

The  diffusion  velocities  {V j }  in  Eqs.  (17,  18)  are', 
determined  by  solving  the  equations  [1,2] 

M  n  .n  VP 

sj  -A  -wr<2k-ij)  =  i  <vN)  -  (p/p  - n/N)^  C19) 

k  — jl  Jr* 

k*  j 

subject  to  the  constraints 

H  H 

T  S.  =  0,  and  l  p.V.  =  0.  (20) 

J-l  J  j=l  3  3 

An  iterative  algorithm  for  solving  for  the  diffusion  veloci¬ 
ties  which  avoids  the  cost  of  performing  matrix  inversions  has 
been  developed  [2,3].  This  algorithm  DPLUX,  is  of  0(M2)  and 
is  vectorized.  Thus  it  is  substantially  faster  than  0(M3) 
matrix  inversions  when  four  or  more  species  are  involved- 

Equations  for  the  evaluation  of  diffusive  transport  coef¬ 
ficients  have  been  discussed  in  detail  by  Picone  [7]  and  sum¬ 
marized  in  Oran  and  Boris  [2].  The  forms  given  are  for  mix¬ 
tures  of  neutral  gases.  Their  derivations  are  based  on  or  are 
advancements  of  the  fundamental  work  of  Chapman  and  Cowling 
[8]  and  Hirschfelder,  Curtis  and  Bird  [9].  Representations  of 
the  coefficients  which  are  easily  used  in  reactive  flow  models 
and  not  requiring  the  expensive  inversion  of  matrices  are 
chosen. 
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2.5  Time step  Sp lilting 

In  the  asymptotic  tiir.estep  split  approach  [2],  the  indi¬ 
vidual  terms  in  the  Eqs.  (1-4)  are  solved  independently  as 
described  above  and  then  asymptotically  coupled  together. 

Since  both  chemistry  and  sound  waves  are  usually  stiff  in 
deflagration  problems,  special  care  is  required  in  coupling 
the  chemical  heat  release  to  the  fluid  dynamics.  In  a  flame, 
fluid  dynamic  expansion  and  diffusive  transport  relieve  the 
pressure  from  the  flame  region  as  fast  as  it  is  generated. 

Thus  the  pressure  stays  effectively  constant.  Small  pressure 
fluctuations,  0(v2/c2),  do  exist  and  are  just  large  enough  to 
drive  the  flows  which  reapportion  the  energy  released  by 
chemistry  or  transported  by  diffusion.  The  chemistry  step 
should  be  taken  at  constant  pressure,  but  it  may  also  be  taken 
at  constant  volume  with  temperature  held  fixed  if  the  profiles 
only  change  slowly  per  timestep.  At  the  completion  of  the 
chemistry  integration,  the  heat  release  is  converted  to  an 
effective  pressure  change  at  constant  volume  because  the  cell 
volume  has  been  held  fixed  during  the  chemical  kinetics 
calculation.  This  pressure  change  is  then  used  as  an  energy 
source  term  for  ~he  fluid  dynamics  timestep  to  get  the  correct 
fluid  expansion.  This  is  done  by  modifying  the  cell  entropies 
used  as  the  fluid  dynamics  module,  ADINC,  described  earlier. 

At  this  point  thermal  conduction  and  diffusion  heat  fluxes 
also  contribute  changes  to  the  entropy. 

2.6  Open  Boundary  Condition 

In  the  study  of  unconfined  flames,  a  representation  for 
the  open-boundary  is  required  since  the  size  of  the  computa¬ 
tional  domain  is  limited.  One  approach  is  to  allow  the  compu¬ 
tational  cells  to  increase  in  size  as  they  get  further  from 
the  flame.  Thus  the  computational  domain  is  very  large  and 
there  is  no  corresponding  increase  in  computer  storage.  How¬ 
ever,  the  cell  stretching  should  be  gradual  to  limit  inaccu¬ 
racies  which  arise  as  a  result  of  the  varying  cell  size. 

A  different  approach  used  in  the  study  of  the  propagation 
of  unconfined  flames  is  discussed  below.  Here  the  effects  of 
an  open  boundary  are  accounted  for  by  allowing  the  pressure 
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changes  which  occur  in  each  cell  during  the  timestep  to  relax 
adiabatically  to  the  pressure  before  the  time  step.  Due  to 
this  relaxation,  the  volume  of  the  cell  changes  since  a 
Lagrangian  coordinate  system  has  been  used.  Tne  changes  in 
the  volume  of  the  cells  causes  the  location  of  the  open  bound¬ 
ary  to  change.  The  location  of  the  open  boundary  and  the 
fluid  velocity  in  the  last  cell  (which  is  also  the  velocity  of 
the  o  en  boundary)  are  used  as  open  boundary  conditions.  This 
procedure  is  appropriate  because  in  unconfined  flames  the 
pressure  stays  effectively  constant. 

3.  FLAMES  IN  HYDROGEN-OXYGEN-NITROGEN  MIXTURES 

Simulations  of  flames  in  hydrogen-oxygen-nitrogen  mix¬ 
tures  have  been  carried  out  using  the  numerical  model  de¬ 
scribed  above.  Applying  the  model  to  a  specific  gas  mixture 
requires  knowledge  of  the  chemical  kinetic  rate  scheme  and 
other  species  dependent  input  parameters  such  as  thermal  con- 
ducitvity  and  molecular  diffusion  coefficients. 

The  chemical  kinetics  rate  scheme  used  consists  of  about 
fifty  chemical  rates  relating  the  species  H2,  02>  H,  0,  OH, 
H20,  H02  and  H2C..  It  has  been  extensively  tested  against 
experimental  data  [10,  113  and  shown  to  give  good  results. 
Burks  and  Oran  [1C  I  showed  that  the  results  computed  with  the 
scheme  compared  very  well  with  experimentally  observed  induc¬ 
tion  times,  second  explosion  limits  and  the  temporal  behavior 
of  reactive  species.  Oran  et  al.  [11]  have  shown  that  the 
scheme  gives  good  results  when  coupled  with  a  fluid  dynamic 
model  in  the  'simulation  of  the  conditions  behind  a  reflected 
shock.  The  reaction  rate  scheme  has  not  been  presented  here 
since  it  is  readily  available  [10,11,123.  Heats  of  formation 
and  enthalpies  have  been  taken  from  the  JANAF  tables  [13]. 

The  diffusion  coefficients  have  been  obtained  using  the 
data  given  by  Marrero  and  Mason  [1*1].  The  Lennard-Jones 
(12:6)  potential  parameters  c  and  ( e/lt)  are  those  given  by 
Svohla  [153  except  for  ojj  and  oJl20  for  which  the  values 
given  by  Dixon -Lewis  [16]  have  boon  used. 

Below  we  first  discuss  some  of  the  results  obtained  in  a 
study  of  the  ignition  and  quenching  of  a  2:1 :10/II2:0  2:N2 
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mixture.  Then  we  present  results  from  a  simulation  of  a 
propagating  flame  in  a  H2-air  mixture.  In  both  cases,  the 
initial  temperatures  and  pressures  of  the  unburned  gas  were 
298  K  and  1  atm,  respectively. 

3.1  Minimum  Ignition  Energies  and  Quenching  Distances 

The  model  described  above  is  ideally  suited  for  studying 
time-dependent  problems  such  as  the  ignition  of  gas  mixtures 
[12].  The  model  was  configured  in  spherical  geometry  with  an 
open  boundary  at  one  end  to  simulate  a  very  large  system. 
Energy  was  deposited  in  the  center  with  a  radius  of  deposi¬ 
tion,  Rq.  Results  from  a  typical  calculation  are  presented 
in  Fig.  1.  The  figure  depicts  the  time  history  of  the  temper¬ 
ature  profile  after  4  mJ  of  energy  is  deposited  over  a  period 
of  0.1  ms.  Even  after  the  energy  deposition  is  stopped,  the 
central  temperature  continues  to  increase  due  to  the  heat 
released  in  chemical  reactions.  With  time,  however,  the  tem¬ 
perature  near  the  center  decreases  and  the  temperature  away 
from  the  center  increases  due  to  diffusive  transport.  By 
4.5  ms,  we  see  that  the  temperature  distribution  exhibits  a 
"flame"  temperature  profile.  If  the  amount  of  energy  depos¬ 
ited,  E0,  was  reduced  to  3  mJ,  the  temperature  distribution 
does  not  develop  into  a  flame  temperature  profile. 


F-ig.J.  Time  hei  toKij  0  the  tenpefi- 
atuAc  pioii.Ce. 


Fig.  2.  liiniwuin  ignition  cue Agtj  <x&  a 
jjuucf.teii  fiadiui  o&  zneAgj 
rlopo  Action. 
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By  repeating  the  computations  for  different  values  of 
Ec,  a  bound  for  the  minimum  ignition  energy  for  that  partic¬ 
ular  radius  was  obtained.  Similar  calculations  were  per¬ 
formed  for  different  values  of  the  radius  of  deposition,  RQ* 

The  results  of  such  investigations  are  shown  in  Fig  2.  A 
propagating  flame  results  when  3*8  mJ  of  energy  is  deposited 
in  a  sphere  with  a  radius  of  0.1  cm.  However  if  the  same 
amount  of  energy  is  deposited  in  a  sphere  of  smaller  radius, 
the  rate  of  heat  liberation  is  insufficient  to  compensate  for 
the  rate  of  heat  loss  and  consequently  there  is  no  ignition. 
This  radius,  0.1  cm,  is  the  "quench-radius"  for  this  partic¬ 
ular  mixture.  For  radii  slightly  larger  than  the  quench- 
radius,  the  minimum  ignition  energy  is  almost  constant  and  for 
larger  radii  (larger  than  0.11  cm)  the  minimum  ignition  energy 
increases  rapidly  with  increasing  radii.  Therefore  for  the 
system  under  study,  the  absolute  minimum  energy  is  about  3*7 
mJ.  These  observations  are  in  qualitative  agreement  with 
those  of  Lewis  and  von  Elbe  [17]*  Quantitative  comparisons 
are  not  possible  since  the  composition  of  the  mixture  and  the 
time  for  energy  deposition  are  different. 

3.2  Calculation  of  the  Burning  Velocity  of  an  H.,:0.>:N-,/2:-l :  l\ 
Mixture  at  29SK,  1  atm. 

The  model  described  above  can  also  be  used  to  study  the 
propagation  of  flames  in  pre-mixed  gas  mixtures.  As  shown  in 
Fig.  1,  if  sufficient  energy  is  deposited  at  the  center  of  a 
sphere  and  computations  are  carried  out  for  a  sufficiently 
long  time,  the  temperature  distribution  attains  a  typical 
"flame  profile”.  However,  this  method  is  an  expensive  way  to 
generate  a  propagating  flame  since  so  much  time  is  spent  in 
establishing  and  overcoming  the  Initial  conditions.  Another 
method  for  initializing  the  problem  quickly  is  to  start  the 
computations  with  a  good  guess  for  the  temperature  and  species 
profile  in  a  region  behind  a  flame  front.  The  closer  the 
initial  profiles  are  to  the  final,  steady  state  flame  profile, 
the  sooner  the  initial  conditions  relax  to  the  steady  propa¬ 
gating  flame.  This  procedure  was  adopted  for  obtaining  the 
values  for  a  flame  velocity. 
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The  model  was  set  up  in  cartesian  Geometry  with  one  end 
closed  (preventing  any  gas  flow  through  it)  and  the  other  open 
to  the  atmosphere.  At  the  closed  end,  a  first  approximation 
to  the  temperature  and  major  species  profile  was  set  up  as 
initial  conditions.  This  soon  evolves  into  a  propagating 
flame  as  can  be  seen  in  Fig.  3>  where  the  temperature  profiles 
at  50  ps  and  60  ps  (from  the  start  of  the  computations  )  are 
shown.  However,  in  order  to  determine  the  flame  speed,  a 
criterion  for  the  location  of  the  flame  is  required.  An  arbi¬ 
trary  value  of  900K  is  chosen  to  define  the  location  of  the 
flame.  The  movement  of  the  location  of  this  value  (900K)  in 
time  and  space  gives  the  rate  of  propagation  of  the  flame. 

The  location  of  the  flame  is  presented  as  a  function  of  time.— 
in  curve  (a)  in  Fig.  M.  The  slope  of  this  curve  gives  the 
flame  velocity  (for  an  inertial  observer).  From  this  figure 
we  see  that  the  flame  propagates  initially  at  a  velocity  of  ' 
8.7  m/s  and  by  50  ps  (from  the  start  of  the  computations)  it 
attains  a  nearly  constant  value  of  9-7  m/s. 


Fig .  3 .  TcwipcAatcAC  pro£iZe&  in  0. 
propagating  *laine. 


Tig. 4.  Time,  klitomj  0$  the.  location 
ojj  (a)  the.  i (Came  and  lb)  the 
open  boundary. 


A  quantity  of  considerable  practical  interest  is  the 
burning  velocity  which  is  defined  as  the  rate  at  which  the 
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flame  consumes  the  reactants,  or  equally  well  as  the  flame 
velocity  relative  to  the  unburned  gases.  In  Fig.  5  the  fluid 
velocity  profile  in  the  system  is  depicted  at  a  particular 
time.  The  velocity  is  zero  at  one  end  since  that  end  is  a 
closed  boundary.  The  flow  velocity  increases  across  the  flame 
and  is  maximum  at  the  open  boundary.  An  estimate  for  the 
burning  velocity  can  be  obtained  by  subtracting  the  velocity 
of  the  unburned  gases  from  the  flame  velocity.  An  upper  esti¬ 
mate  for  the  velocity  of  the  unburned  gases  is  the  fluid 
velocity  at  the  open  boundary.  Since  this  model  uses  a 
Lagrangian  coordinate  system,  the  fluid  velocity  at  the  open 
boundary  is  also  given  by  the  rate  of  movement  cf  the  open 
boundary.  The  location  of  the  open  boundary,  has  been  depicted 
as  a  function  of  time  in  Fig.  4b.  Corresponding  to  a  flame 
velocity  of  9-7  m/s  we  see  that  the  velocity  of  the  unburned 
gases  (open-boundary)  is  about  7.6  m/s  giving  a  burning  veloc¬ 
ity  of  2.1  m/s  for  the  hydrogen-air  flame  studied.  This  com¬ 
pares  well  with  the .  experimental  values  for  the  burning  veloc¬ 
ity  Cl 8 3  which  are  between  2  and  2.5  m/s. 


Fig.  S.  Flow  velocity  pxofcZe  in  a 
propagating  Itamc. 


IV.  CONCLUSIONS 

In  this  paper  we  have  shown  that  the  time-dependent  model 
described  does  a  good  job  in  predicting  both  transient  and 


steady  state  pruperti.es  of  laminar  hydroEcn-oxygen-nitrogen 
flames.  In  the  calculations  of  the  burning  velocity  of  a 
mixture  of  H2 :02  :N2/2: 1 :  *1  at  1  atm  and  298  K,  the  velocity 
predicted  was  at  the  lower  end  of  the  observed  range.  We  note, 
however,  that  there  are  several  factors  which  would  tend  to 
increase  this  estimated  value.  First,  the  calculated  burning 
velocity  would  be  larger-  if  the  velocity  of  the  unburned  gas 
used  in  this  calculation  were  closer  to  the  flame  front  in¬ 
stead  of  aL  the  open  boundary.  It  is  not  clear  where  it 
should  be  evaluated.  Second,  if  thermal  diffusion  is  in¬ 
cluded,  it  would  increase  the  flame  velocity  by  perhaps  5-10JS. 
Finally,  the  extremely  important  hydrogen  diffusion  coeffi¬ 
cients  are  not  really  known  at  high  temperatures.  If  they  are 
higher  than  the  ones  used  here,  the  flame  velocity  would 
increase.  These  observations  lead  us  to  believe  that  the 
higher  observed  values  of  the  burning  velocity  are  probably  ■ 
more  accurate. 

Finally,  we  wish  to  note  that  the  method  is  sensitive  to 
initial  conditions.  As  in  the  steady-state  calculations,  a 
better  initial  guess  leads  to  faster  convergence  to  a  steady 
propagating  flame.  Perhaps  the  best  way  to  initialize  such  a 
problem  is  to  perform  a  two  step  process  in -which  we  first 
determine  the  minimum  ignition  energy  for  a  given  radius  and 
time  of  energy  deposition,  and  then  deposit  this  energy  and 
evaluate  the  properties  of  the  propagating  flamefront.  How¬ 
ever,  considering  that  the  experimental  velocities  vary  by 
~252,  a  good  Initial  guess  at  temperature  and  species  concen¬ 
trations  certainly  gives  an  adequate  calculation.  Finally,  we 
note  that  the  criterion  on  which  to  judge  the  initial  guess 
values  is  to  look  at  how  flat  the  temperature  profile  is 
behind  the  flamefront.  A  bad  guess  produces  some  structure 
and  oscillations  which  only  slowly  disappear. 
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ABSTRACT 


Detailed  numerical  simulations  of  shock  tube  exper: 
ments  have  been  used  to  study  the  direct  initiation  of 
detonations  in  hydrogen-air  mixtures.  Whereas  these 
simulations  qualitatively  agree  Kith  earlier  experimental 
observations  that  a  minimum  power,  a  minimum  energy  and 
some  relation  between  them  is  required  to  characterize  the 
direct  initiation  of  detonations,  the  experiments  and 
calculations  do  not  agree  well  quantitatively.  During  the 
course  of  this  study  it  was  found  that  the  disagreement  is 
primarily  because  the  induction  times  of  the  gaseous  mix¬ 
tures,  at  the  temperatures  and  pressures  expected  behind 
the  incident  shocks,  are  very  sensitive  to  perturbations 
in  pressure  and  temperature.  Simulations  carried  out  to 
evaluate  the  effects  of  sound  wave  perturbations  on  one  of 
the  gas  mixtures  studied  give  a  quantitative  relation 

» 

between  the  induction  time  and  sound  wave  perturbations 
in  a  range  of  amplitudes  and  frequencies.  It  is  shown 
that  the  presence  of  such  perturbations  can  produce  experi 

» 

mentally  determined  minimum  powers  and  minimum  energies 
which  are  significantly  different  from  those  calculated. 


I.  INTRODUCTION 


The  widespread  use  of  natural  gas,  methane  and 
other  gaseous  fuels  has  made  it  increasingly  important  to 
understand  and  quantify  their  tendency  to  detonate.  The 
early  works1,2,5  on  the  initiation  of  detonations  estab¬ 
lished  the  importance  of  the  magnitude  of  the  source 
energy  in  the  direct  initiation  of  gaseous  detonations. 

More  recent  experiments4,5’^1  have  shown  that  not  only  is 
the  amount  of  energy  important  but  also  the  rate  at  which 
the  energy  is  deposited,  namely,  the  power.  The  experi¬ 
mental  results  cf  Lee  et  al.5  indicate  that  there  is  a 
critical  detonation  energy  below  which  a  detonation  would 
not  occur  no  matter  what  the  power  is  and  that  there 
is  a  critical  power  below  which  a  detonation  would  not 
occur  no  matter  what  the  total  energy  is.  The  require¬ 
ment  for  a  critical  value  for  the  power  of  the  source 
indicates  that  the  source  must  be  capable  of  generating  a 
shock  wave  of  certain  minimum  strength.  The  critical  energy 
requirement  implies  that  the  shock  wave  must  be  maintained 
at  or  above  this  minimum  strength  for  a  certain  duration. 

Recently  these  ideas  have  been  used  to  obtain  a 
relation  between  the  power  and  energy  requirement  for  direct 
initiation  of  detonations  in  a  shock  tube.'  In  a  shock 
tube,  for  a  given  set  of  driver  and  driven  gases,  the 
pressure  ratio  across  the  diaphragm  determines  the  strength 
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(Mach  Number)  of  the  shock  wave  generated.  Therefore,  by 
varying  the  pressure  ratio  across  the  diaphragm,  the  power 
of  the  source  can  be  varied  since  it  depends  on  the  shock 
strength.  When  the  diaphragm  is  burst,  rarefaction  waves 
are  generated  in  the  driver  section  of  the  shock  tube.  If 
a  short  driver  section  is  used  the  rarefaction  waves  can 
reach  the  shock  front  after  reflection  from  the  driver  end 
wall  and  reduce  the  shock  velocity.  If  the  wave  is  a 
detonation,  it  is  expected  that  the  wave  will  continue  to 
propagate  unhampered  by  the  rarefaction  waves.  For  a  specific 
driver  length,  we  can  determine  the  pressure  ratio  (and 
hence  the  minimum  power)  below  which  the  rarefaction  waves 
will  interact  with  the  shock  wave  and  cause  it  to  decay 
before  ignition  can  occur.  The  corresponding  minimum 
detonation  energy  can  be  obtained  if  the  time  for  which 
the  shock  wave  is  effective  is  known.  By  repeating  the 
experiments  with  different  driver  lengths,  we  can  determine 
different  sets  of  minimum  power  and  minimum  energy  required 
to  initate  detonations.  Then,  from  a  plot  of  these  minimum 
powers  as  a  function  of  the  minimum  energies,  one  can  deter¬ 
mine  the  critical  power  (the  lowest  value  of  the  minimum 
power)  and  the  critical  detonation  energy  (the  lowest 
value  of  the  minimum  energy). 
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In  the  shock  tube  experiments,  the  initiation  of 
stoichiometric  hydrogen-air  mixtures  with  drivers  of  three 
different  lengths  was  investigated.  For  each  driver,  the 


location  at  which  the  Chapman- Jouguet  Mach  Number  was 
reached  was  noted  for  various  pressure  ratios.  From  this 
the  pressure  ratio  below  which  detonation  was  expected  to 
take  place  at  extremely  long  distances  from  the  diaphragm 
was  estimated.  Using  this  information  and  a  simple 
theoretical  model  a  relation  between  the  minimum  detonation 
energy  and  the  minimum  power  was  also  obtained. 

In  this  paper,  we  have  used  numerical  simulations 
to  examine  several  aspects  of  the  basic  idea  of  using  shock 
tube  experiments  to  study  the  direct  initiation  of  detona¬ 
tions.  First,  we  establish  the  importance  of  the  power  in 
the  direct  initiation  of  detonations.  Then  we  show  that 
the  results  from  the  numerical  simulations  are  very  dif¬ 
ferent  from  the  experimental  results.  The  reasons  for 
the  discrepancy  are  examined  and  discussed  in  detail. 

First,  the  experiments  are  shown  to  be  in  a  pressure- 
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temperature  regime  which  is  very  sensitive  to  perturbations  * 
and  second,  the  experiments  were  very  susceptible  to  such 
perturbations  because  the  driver  and  driven  sections  had 
different  areas  of  cross-section  and  were  of  different 
shapes.  Finally,  we  examine  this  sensitivity  of  the  igni¬ 
tion  process  in  the  hydrogen-air  mixture  by  using  numerical 
simulations  to  evaluate  the  effect  of  sound  wave  perturbations 
on  ignition. 
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II.  THE  NUMERICAL  MODEL 

The  one-dimensional  N'RL  Reactive  Shock  Model1®’** 
used  to  perform  the  calculations  described  below  solves 
the  time -dependent  conservation  equations** for  mass, 
momentum  and  energy  coupled  to  the  equations  describing  the 
chemical  kinetics.  Though  the  geometry  may  be  either 
cartesian,  cylindrical,  spherical  or  some  generalized  co¬ 
ordinate,  here  we  consider  cartesian  geometry.  The  model 
uses  an  explicit,  Eulerian  finite  difference  formulation 
with  a  sliding  recone  capability  to  provide  resolution 
around  moving  gradients. 

The  convective  transport  terms  in  the  conservation 
equations  are  solved  using  one  variant  of  the  Flux-Corrected 
Transport  (FCT)  method.13’*4  This  is  a  conservative,  mono¬ 
tonic  algorithm  with  fourth-order  phase  accuracy  and  does  not 
require  artificial  viscosity  to  stabilize  shocks.  The 
Flux  Correction  procedure  itself  ensures  that  the  shocks 
are  one  or  two  zones  wide  and  have  maximal  resolution. 

The  ordinary  differential  equations  describing  the  chemical 
kinetics  are  solved  using  YSAIM,  a  vectorized  version  of 
the  selected  asymptotic  integration  method  employed  in 
CHEMEQ.*^’*®  This  algorithm  identifies  the  stiff  equations 
for  treatment  with  a  stiffly  stable  method.  The  remaining 
equations  are  solved  with  a  standard  classical  method. 

The  algorithm  has  been  specially  optimized  for  use  in 
conjunction  with  fluid  dynamic  models. 
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For  calculations  performed  in  this  paper  the  time- 
scales  under  consideration  are  short  and  therefore  the 
diffusive  transport  processes,  thermal  conduction  and 
molecular  di f fusion ,  have  negligible  effect.  The  calcula¬ 
tions  performed  in  this  paper  do  not  include  these  processes 
although  they  are  part  of  the  NRL  Reactive  Shock  Model. 

The  solutions  for  the  equations  describing  the  fluid 
dynamics  and  the  chemistry  of  the  problem  are  coupled 
using  time-step  splitting  techniques.11 

The  chemical  kinetics  rate  scheme  used  consists  of 
about  fifty  chemical  rates  relating  the  species  ,  02, 

H,  0,  OH,  r!C,,  K.O  andH702*  It  has  been  extensively  tested 

O  17 

against  experimental  data  ’  and  shown  to  give  good 
results.  Burks  and  Oran17  showed  that  the  results  com¬ 
puted  with  the  scheme  compared  very  well  with  experimentally 

observed  induction  times,  second  explosion  limits  and  the 

0 

temporal  behavior  of  reactive  species.  Oran  et  al.  have 

shown  that  the  scheme  gives  good  results  when  coupled  with 

a  fluid  dynamic  model  in  the  simulation  of  the  conditions 

behind  a  reflected  shock.  The  reaction  rate  scheme  has  not 
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been  presented  here  since  it  is  readily  available.  ' 

Heats  of  formation  and  enthalpies  have  been  taken  from  the 
JAXAF  tables.18 

The  detailed  simulations  of  the  shock  tube  experi¬ 
ments  discussed  in  this  paper  require  that  we  not  only 


1 
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model  very  long  systems  but  that  we  also  accurately  model 
the  movement  of  the  contact  surface  and  the  reactive  shock 
wave.  A  rather  sophisticated  adaptive  gridding  method  has 

1  Q 

been  developed  for  this  purpose  at  NRL  by  Oran  and  Young. 
There  are  two  finely  gridded  regions:  one  surrounding  and 
moving  with  the  shock  wave  and  the  other  surrounding  and 
moving  with  the  contact  surface.  Each  of  these  finely 
gridded  regions  may  have  a  different  minimum  computational 
cell  sice.  The  regions  ahead  of  the  shock  wave  and  behind 
the  contact  surface  have  computational  grids  in  which  the 
cells  change  exponentially  in  size  from  the  smallest  sur¬ 
rounding  the  shock  wave  or  the  contact  surface  to  the  largest 
at  the  walls.  Care  is  taken  that  the  transition  in  the 
cell  sites  is  smooth.  For  the  results  presented  in  this 
paper  a  total  of  165  computational  cells  are  used  to 
describe  the  shock  tube  and  the  cell  sites  varied  from 
0.1  cm  around  the  shock  to  over  50  cm  near  the  shock  tube 


end-walls. 


III.  NUMERICAL  SIMULATIONS 


The  numerical  model  described  in  the  preceeding 

section  was  first  used  to  study  a  selected  group  of  the 

shock  tube  experiments  described  briefly  in  the  introduc- 

7 

tion  and  more  completely  by  Dabora.  In  the  experiments, 
the  driven  section  of  the  shock  tube  was  9  feet  (“2.75m) 
long  with  2"  x  2"  square  cross  section.  The  driven  section 
contained  a  hydrogen-air  mixture  at  0.5  atm  (50.66  kPa) . 

The  driver  section  contained  helium,  was  variable  in  length 
(2”,  4"  or  8")  and  had  a  circular  cross  section  providing  a 
driver  to  driven  area  ratio  of  1.65. 

The  simulations  described  below  are  for  a  shock 
tube  with  a  20  cm  (»8")  driver  section  filled  with  helium 
and  a  280  cm  driven  section  filled  with  a  mixture  of 
H2-02-N2/2- 1-4  at  0.5  atm  and  298  K.  However,  in  contrast 
to  the  experimental  set-up,  the  entire  tube  was  assumed  to 
be  of  uniform  cross  section. 

One  of  the  problems  with  numerical  simulations  is 
interpreting  the  large  quantity  of  output  data.  At  each 
time  step  we  calculate  mass,  momentum,  energy,  temperature, 
and  a  number  of  species  densities.  In  order  to  compare  the 
simulations  most  effectively  with  the  experiment,  we  have 
chosen  to  follow  the  time-dependent  behavior  of  three 
quantities:  the  velocity  of  the  shock  or  detonation,  the 

velocity  of  the  contact  surface  and  the  pressure  behind 


the  leading  edge  of  the  shock  or  detonation  wave.  Typical 
values  of  these  quantities  are  shown  in  Fig.  1  for  a  pres¬ 
sure  ratio  of  110  across  the  diaphragm.  After  the  diaphragm 
is  burst,  the  velocity  of  the  contact  surface  (Fig-  lb)  is 
almost  constant  until  about  0.6  ms  and  then  it  decreases 
continuously  because  of  interactions  with  the  rarefaction 
waves  reflected  from  the  driver  end-wall.  When  the  front 
of  the  rarefaction  waves  finally  catches  up  with  the  shock 
front,  the  shock  velocity  and  pressure  are  reduced.  The 
shock  does  not  transition  to  detonation  during  the  time  (1.6  ms) 
for  which  the  simulation  was  carried  out. 

Results  from  a  simulation  corresponding  to  a  pres¬ 
sure  ratio  of  150  are  depicted  in  Fig.  2.  Here  ignition 
was  first  observed  at  about  0.14  ms  near  the  contact  sur¬ 
face.  The  chemical  energy  released  causes  a  rapid  rise  in 
the  temperature  and  the  pressure  near  the  contact  surface 
whose  velocity  then  changes  rapidly.  The  entire  region 
between  the  contact  surface  and  the  shock  front  had 
ignited  by  0.15  ms.  The  pressure  behind  the  shock  front 
increases  greatly  at  about  this  time  and  then  decreases 
towards  a  constant  value.  The  shock  transitions  to  a 
detonation  around  0.15  ms.  There  are  oscillations  in  the 
shock  velocity  during  the  transition  period  as  well  as  for 
a  short  time  after  the  diaphragm  is  burst.  However,  by 
about  0.2  ms  the  oscillations  die  down  and  the  detonation 
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velocity  gradually  decreases  towards  the  Chapman- Jouguet 
detonation  velocity.  During  the  time  oscillations  were 
observed,  a  least  squares  approximation  was  used  to  obtain 
a  fit  for  the  shock  or  detonation  front  velocity. 

These  calculations  can  be  used  to  evaluate  the 
power,  which  is  the  rate  of  energy  input  to  the  driven 
section.  The  power  is  equal  to  the  rate  of  work  performed 
by  the  contact  surface  which  is  the  product  of  the  pressure 
at  the  contact  surface  and  the  velocity  of  the  contact 
surface.  The  pressure  at  the  contact  surface  is  directly 
obtained  from  the  simulations  while  the  velocity  of  the 
contact  surface  is  calculated  from  the  time  history  of  the 
location  of  the  contact  surface. 

The  time  history  of  the  power  and  energy  for  the 
two  cases  discussed  in  detail  above  are  shown  in  Fig.  3. 

In  the  case  where  the  pressure  ratio  was  130,  the  power 
is  nearly  constant  until  about  0.6  ms  and  then  it  decreases 
as  the  pressure  and  the  fluid  velocity  at  the  contact 

surface  decrease.  Figure  3a  shows  that  the  power  decreases 

2  2 

from  about  0.9  GK/rn  to  about  0.34  GIV/m  by  1.6  ms.  The 
energy  input  to  the  driven  gas  increases  linearly  until 
0.6  ms  since  the  power  is  constant.  It  then  increases  less 
rapidly  as  the  power  decreases  with  time.  Ke  note  that  an 

7 

energy  of  1.12  MJ/m  is  not  sufficient  to  cause  a  detona¬ 
tion  in  this  case.  In  the  case  where  the  pressure  ratio 
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was  150  (Fig.  5b),  the  power  is  almost  constant  until  0.14  ms, 

rises  rapidly  and  fluctuates  until  about  0.24  ms,  and  then 

remains  constant  until  the  end  of  the  simulation.  The  rapid 

rise  in  the  power  is  due  to  the  interaction  of  the  contact 

surface  with  pressure  waves  generated  by  combustion.  The 

significance  of  the  power  after  ignition  occurs  is  not 

clear  since  then  the  energy  released  in  chemical  reactions 

interact  with  the  motion  of  the  contact  surface.  The  energy 

2 

input  to  the  driven  gases  up  to  ignition  is  only  0.135  MJ/m  . 

In  summary,  we  observe  from  these  simulations  that 
the  0.9  GV.'/m"  initial  power  with  an  energy  input  of 
1.12  MJ/m"  does  not  initiate  a  detonation  while  the 

9  2 

1.1  Gh'/m"  initial  power  corresponding  to  0.135  MJ/m  energy 
input  does.  This  indicates  that  both  a  minimum  power, 
a  minimum  energy  and  some  relation  between  them  are 
required  to  characterize  the  initiation  of  detonations. 
However,  before  any  further  effort  goes  into  determining 
them,  it  is  important  to  validate  this  approach  by  comparing 
the  calculated  and  experimental  results. 

The  pressure  ratios  and  the  corresponding  detona¬ 
tion  locations  observed  in  the  experiments  and  in  the 
simulations  are  shown  in  Fig.  4,  Both  the  simulations  and 
the  experiments  show  that  as  the  pressure  ratio  is  decreased, 
detonation  occurs  further  and  further  from  the  diaphragm. 
However,  the  numerically  predicted  pressure  ratio  below 
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which  a  detonation  will  occur  far  from  the  diaphragm  as 
well  as  the  locations  of  the  detonations  are  very  different 
from  the  experimental  results.  Part  of  the  difference 
between  the  experimentally  and  numerically  determined 
locations  is  due  to  different  criteria  used  to  determine 
them.  The  experimental  locations  are  determined  by  the 
observation  of  a  "steady,"  Chapman-Jouguet  detonation 
velocity  while  the  numerically  determined  locations  are 
those  at  which  an  "unsteady"  detonation  wave  is  first 
observed.  In  the  simulations,  the  detonation  velocity  was 
found  to  be  gradually  decreasing  towards  a  steady  value 
when  the  computations  were  terminated.  Longer  runs  are 
required  to  determine  the  location  at  which  the  detonation 
wave  becomes  steady. 

A  possible  explanation  for  the  differences  in  the 
pressure  ratio  for  which  detonations  are  observed  could  be 
the  difference  in  the  driver  to  driven  area  ratio.  For  a 
specific  pressure  ratio,  the  higher  area  ratio  in  the 
experiments  may  result  in  a  shock  of  higher  Mach  number 
compared  to  those  simulated.  However,  the  Mach  numbers  of 
the  shocks  actually  observed  in  the  experiments  (Fig.  4) 
were  all  lower  than  those  in  the  simulations.  Therefore 
an  alternate  explanation  is  required  for  the  discrepancy 
between  the  two  results. 


The  disagreement  between  the  predictions  of  the 

numerical  model  and  the  experimental  results  is  similar  to 

that  observed  in  an  earlier  study  of  weak  ignition  systems 
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by  Oran  and  Boris.  *  They  found  that  ignition  in  mixtures 
of  hydrogen  and  oxygen  in  certain  regimes  of  the  pressure- 
temperature  plane  are  extremely  sensitive  to  sound  wave  and 
entropy  wave  perturbations.  In  the  next  section,  the 
discrepancy  in  the  results  is  examined  on  the  basis  of 
this  observation. 


IV.  SENSITIVITY  OF  IGNITION  TO  PERTURBATIONS 


The  sensitivity  of  ignition  to  perturbations  can  be 
assessed  by  considering  derivatives  of  the  chemical  induction 
time  vith  respect  to  parameters  of  the  system  such  as  temper¬ 
ature,  pressure,  or  stoichiometry.  One  important  quantity 
to  consider  is 

S  =  - (T/t) (3t/3T I $) > 

which  shows  the  sensitivity  of  the  induction  time,  t ,  to 

9 

sound  wave  perturbations.  Here  T  and  s  are  the  initial 

q 

temperature  and  entropy,  respectively.  It  has  been  shown' 
that  high  values  of  S  indicate  large  variations  in  the 
ignition  time  for  small  changes  in  pressure  and  temperature. 

The  value  of  S  for  the  conditions  behind  the  incident  shock 
in  the  calculations  discussed  above  were  all  between  35  and 
40,  indicating  a  fairly  strong  susceptibility  of  ignition 
times  to  perturbations.  This  means  that  2.5-3%  change  in 
temperatures  changes  the  induction  time  by  100%. 

In  order  to  evaluate  quantitatively  how  a  specific 
type  of  perturbation  affects  ignition,  we  have  simulated 
the  effects  of  sound  waves  in  H2-0,-N2  mixtures  by  recon¬ 
figuring  the  numerical  model  described  earlier.  The  particular 
case  we  will  investigate  here  is  that  for  the  pressure  ratio 
of  110,  which  corresponds  to  a  high  value  of  S  in  the  cal¬ 
culations  considered  in  the  previous  section,  h'e  consider 
a  homogeneous  mixture  of  H,-02'N->/-- 1- 4  at  975K  and  800kPa. 
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At  the  beginning  of  these  simulations,  a  velocity  perturbation 
is  imposed  on  the  system  at  each  location  x  such  that 

v(x,t=0)  =  vosin(— j~—  x) 

where  vq  is  the  amplitude  and  L  is  the  half  wavelength  of 
the  sound  wave  perturbation.  We  determine  L  by  deciding  how 
many  periods  of  the  wave  we  wish  to  occur  during  a  chemical 
induction  time. 

Figure  5  shows  the  time  history  of  the  temperature 
taken  from  one  simulation  in  which  there  are  roughly  three 
periods  of  sound  wave  oscillation  in  an  induction  time.  In 
this  case,  L  is  90  cm  (giving  a  frequency  of  405  Hz)  since  the 
induction  time  is  7580}js.  The  amplitude  of  the  perturbation,  vq, 
is  lOOm/s.  With  this  perturbation,  the  system  ignites  first  at 
the  left  wall  at  4925us,  then  at  the  center  and  finally  at  the 
right  wall  at  5091ys.  Thus  we  observe  that  the  time  to  ignition 
is  reduced  by  2457ps  due  to  the  sound  wave  perturbation. 

Figure  6  summarizes  the  results  of  a  number  of 
calculations  similar  to  the  one  presented  in  Fig.  5.  Here 
we  have  plotted  the  time  at  which  ignition  is  first  observed 
as  a  function  of  vQ  for  a  range  of  values  of  L.  For  example, 
if  L  is  90cm,  perturbations  with  amplitudes  greater  than  lOOm/s 
reduce  the  induction  time  substantially.  Considering  the 
curves  for  two  values  of  L,  namely  45cm  and  90cm,  we  note 
that  for  low  amplitude  perturbations  (less  than  250m/s)  ,  the 
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lover  frequency  perturbations  cause  ignition  earlier  than 
the  higher  frequency  perturbation.  The  trend  is  reversed 
for  higher  amplitude  perturbations  Cgreater  than  250m/s) . 

This  difference  can  be  explained  by  considering  the  number 
of  periods  of  oscillation  that  occur  before  ignition.  For 
the  higher  amplitude  perturbations,  the  detailed  simulations 
show  that  ignition  occurs  during  the  first  period  of  the 
perturbations.  Since  the  period  of  oscillation  is  shorter 
with  the  high  frequency  perturbations,  ignition  occurs  earlier 
in  that  case.  For  lover  amplitude  perturbations,  the  system 
undergoes  more  than  one  period  of  oscillation  before  ignition. 
Although  the  period  itself  is  longer  for  lover  frequency 
perturbations,  the  residence  time  for  the  reactive  mixture 
in  a  perturbed  state  is  also  longer.  Therefore,  for  low 
amplitude  perturbations,  when  more  than  one  period  of 
oscillation  occurs  before  ignition,  the  lover  frequency 
perturbation  causes  earlier  ignition.  This  competition 
between  the  residence  time  in  a  perturbed  state  and  the 
period  of  oscillation  explains  the  intersection  of  the 
curves  in  Fig.  6. 

he  also  observe  from  Fig.  6  that  low  frequency 
perturbations  are  more  critical,  since  for  these  frequencies 
smaller  amplitude  perturbations  are  sufficient  to  substantially 
reduce  the  ignition  times.  For  example,  to  reduce  the  time  for 
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ignition  by  a  factor  of  three,  we  need  low  frequency 
perturbations  (L  =  270cm)  of  amplitude  lOOm/s  or  higher 
frequency  perturbations  (L  =  2.7cm)  of  amplitude  330m/s. 

[he  also  note  that  in  the  limit  of  no  frequency  (i.e., 
local  hot  spots)  a  small  temperature  perturbation  could 
have  an  enormous  effect  on  the  ignition  time.) 

In  the  previous  section,  we  noted  differences 
between  the  experimental  results  and  the  numerical  predic¬ 
tions  of  the  initiation  of  detonations.  Specifically,  for 
a  pressure  ratio  of  110,  the  numerical  predictions  did 
not  predict  ignition  while  a  Chapman- Jouguet  detonation 
wave  was  observed  in  the  experiments  for  even  smaller  pres¬ 
sure  ratios.  Here  we  have  shown  that  sound  wave  perturba¬ 
tions  can  cause  large  reductions  in  ignition  times.  For 
example,  a  perturbation  of  amplitude  200m/s  and  half¬ 
wavelength  of  90cm  causes  ignition  as  early  as  660us  which 
is  sufficient  to  explain  the  discrepancy  between  the 
numerical  predictions  and  the  experimental  observations. 
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V.  CONCLUSIONS 
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In  this  paper  detailed  numerical  simulations  of 

incident  shock  tube  experiments  have  been  used  first  to 

S  7 

verify  earlier  experimental  observations  ’  that  a  minimum 
power,  a  minimum  energy  and  some  relation  bet\\-een  them  is 
required  to  characterize  the  direct  initiation  of  detona- 
tions.  However,  in  the  course  of  this  study  it  became 
necessary  to  reconcile  the  quantitative  differences  between 
the  numerical  predictions  and  the  experimental  observations. 
This  in  turn  led  us  to  investigate  the  sensitivity  of  the 
mixtures  in  the  region  behind  the  incident  shock  to  fluctua¬ 
tions  in  pressure  and  fc-mperature  which  might  cause  early 
ignition  in  an  experiment. 

Simulations  were  carried  out  to  evaluate  the  effects 
of  sound  wave  perturbations  which  in  turn  result  in  tempera¬ 
ture  and  pressure  fluctuations.  The  results  of  these  sound 
wave  studies  shows  how  the  ignition  time  is  affected  by 
sound  waves  of  a  given  amplitude  and  frequency.  We  note 
that  low  frequency  perturbations  (~480  Hz  or  half  wavelength 
L  ~  75  cm)  of  large  enough  amplitude  (~200  m/s)  can  reduce 
the  ignition  time  by  an  order  of  magnitude.  We  also  observe 


that  low  frequency  perturbations  are  more  likely  than  high 
frequency  perturbations  to  cause  the  observed  discrepancies 
between  the  experiments  and  the  calculations. 

The  presence  of  temperature  and  pressure  variations 
in  the  experiments  can  be  due  to  factors  such  as  the  non-ideal 
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break-up  oi'  the  diaphragm,  boundary  layers  causing  non¬ 
uniformities  in  the  walls  of  the  tube,  the  sudden  change 
in  shape  and  area  of  cross-section  between  the  driver  and 
driven  sections,  or  heat  loss  to  the  walls.  However, 
perturbations  required  to  substantially  alter  the  ignition 
times  to  cause  such  discrepancies  as  observed  here  may  be 
larger  than  what  might  normally  be  present  in  very  well 
controlled  shock  tube  experiments.  This  needs  to  be 
investigated  before  one  can  confidently  use  shock  tube 
studies  to  determine  the  critical  conditions  necessary  to 
characterise  the  direct  initiation  of  detonations. 

Different  systems  will  have  different  temperature, 
pressure,  sound  wave  or  other  perturbations  and  \<e  have 
shown  here  that  such  perturbations  can  substantially  reduce 
the  ignition  time  and  hence  the  critical  conditions.  It 
appears  that  the  critical  conditions  obtained  from  different 
experimental  configurations  will  be  different  from  each 
other  and  from  practical  systems.  Therefore,  further  work 
needs  to  be  done  to  estimate  the  perturbations  likely  to 
occur  in  various  systems. 
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FIGURE  CAPTIONS 


Figure  1. 


Figure  2. 


Figure  5. 


Figure  4. 


Figure  5. 


Figure  6. 


Time  history  of  (a)  the  shock  velocity,  (b)  the 
contact  surface  velocity,  and  (c)  the  pressure 
behind  the  shock  front  for  the  diaphragm 
pressure  ratio  of  110. 

Time  history  of  (a)  the  shock  or  detonation 
velocity,  (b)  the  contact  surface  velocity,  and 
(c)  the  pressure  behind  the  shock  front  for  the 
diaphragm  pressure  ratio  of  150.  The 
Chapman- Jouguet  detonation  velocity  is  2  km/s. 
Time  history  of  the  power  and  energy  for  the 
diaphragm  pressure  ratio  of  (A)  110  and  (B)  150. 
The  pressure  ratios  and  the  corresponding  detona 
tion  locations  observed  in  the  experiments  and 
in  the  numerical  simulations.  M  is  the  initial 
Mach  number  of  the  shock. 

Time  history  of  the  temperature  at  three  loca¬ 
tions  in  a  0.9m  long  system  with  the  sound  wave 
perturbation  of  amplitude  100  m/s. 

The  ignition  time  as  a  function  of  the  amplitude 
of  the  sound  wave  perturbation  in  systems  of 
various  length. 
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